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The multiple gravity assist low-thrust (MGALT) trajectory model combines the medium-fidelity 
Sims-Flanagan bounded-impulse transcription with a patched-conics flyby model and is an important 
tool for preliminary trajectory design. While this model features fast state propagation via Kepler’s 
equation and provides a pleasingly accurate estimation of the total mass budget for the eventual flight- 
suitable integrated trajectory it does suffer from one major drawback, namely its temporal spacing of 
the control nodes. We introduce a variant of the MGALT transcription that utilizes the generalized 
anomaly from the universal formulation of Kepler’s equation as a decision variable in addition to 
the trajectory phase propagation time.This results in two improvements over the traditional model. 
The first is that the maneuver locations are equally spaced in generalized anomaly about the orbit 
rather than time. The second is that the Kepler propagator now has the generalized anomaly as its 
independent variable instead of time and thus becomes an iteration-free propagation method. The 
new algorithm is outlined, including the impact that this has on the computation of Jacobian entries 
for numerical optimization, and a motivating application problem is presented that illustrates the 
improvements that this model has over the traditional MGALT transcription. 


INTRODUCTION 


Bounded-impulse trajectory models are an integral part of many preliminary design methodologies. For the case of 
continuous-thrust trajectory design, the multiple gravity-assist low-thrust (MGALT) model has been incorporated 
into many tool chains such as the Evolutionary Mission Trajectory Generator (EMTG), Mission Analysis Low- 
Thrust Optimization (MALTO),! Gravity-Assist Low-thrust Local Optimization Program (GALLOP)” and the Parallel 
Global Multiobjective Optimizer (PaGMO).* The MGALT transcription combines the medium-fidelity Sims-Flanagan 
bounded-impulse model,* which divides a trajectory phase into N equal sized time segments with a bounded-impulse 
placed at the center of each segment, with a patched-conics flyby model. Its strength lies in generating accurate mass 
budgets for interplanetary missions. MGALT optimization runs generate medium-fidelity solutions and, therefore, are 
typically used as seeds for higher fidelity models such as the finite-burn low-thrust (FBLT) model®° and flight-grade 
tools such as GMAT’ and Mystic.*-? 

While MGALT is capable of accurately approximating many types of interplanetary trajectories, certain classes of 
missions such as those requiring multiple-revolutions of the central body, with no intermediate flybys, as well as those 
requiring large eccentricity changes present more of a challenge to this model. The more revolutions that the spacecraft 
requires to complete a transfer, the more control discretization nodes are required to ensure that the bounded-impulse 
trajectory will be suitable as a seed to a finite-burn model. For trajectories requiring large eccentricity changes, the 
original Sims-Flanagan model, which evenly spaces control nodes in time, places these nodes at precisely non-optimal 
locations. That is, for high-eccentricity orbits, the control nodes cluster around apoapsis and are sparse at periapsis 
where energy changes are more efficiently achieved and therefore where higher control authority is typically desirable. 
The only way to mitigate this shortcoming of MGALT is to add more control points, which increases the computational 
complexity of the problem. For this reason, a bounded-impulse transcription that requires fewer control nodes and, 
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at a minimum, results in a geometrically even distribution of the control points around the trajectory phase would be 
desirable. 

Redistribution of the control nodes is typically achieved using a time-regularization transformation such as a Sund- 
man transformation. This has recently been investigated by Pellegrini et al.!° who applied generalized Sundman trans- 
formations to the Stark and Kepler (sub-case of Stark) problems. Their Stark solution method is based on a series 
solution, where the length of the series determines the accuracy of the solution. The Stark model was also extended 
for use with a Taylor integrator similar to the work by Yam et al.'' In this paper, we introduce an exact method for 
solving the time-regularized Kepler problem that does not rely on a series solution. Solving Kepler’s equation typically 
requires iterative techniques such as Newton’s method, the Laguerre-Conway method or bisection. We propose a strat- 
egy that transforms a universal Kepler propagator into an iteration-free method at the expense of the introduction of 
one additional decision variable per trajectory phase (to be chosen by an NLP optimizer). These new variables are the 
total generalized anomaly to propagate the spacecraft state by in a trajectory phase (yp). The addition of these decision 
parameters results in an execution speed increase due to the elimination of numerical iteration in the propagator, and 
also a more favorable distribution of control nodes for highly eccentric trajectories. In addition, it is possible to com- 
pute analytic partial derivatives of the problem constraint vector c with respect to the x, variables. The redistribution 
of nodes using this method is a Sundman transformation,!®:!? as such we will hereafter refer to the modified model as 
the multiple gravity-assist low-thrust Sundman (MGALTS) transcription. 


ALGORITHM 


For the case of the original MGALT model, a numerical optimizer selects the phase flight time variable (At,), and a 
Kepler solver then determines the corresponding change in eccentric anomaly, or generalized anomaly, for a universal 
variable formulation. When time-regularization is applied, where control nodes are distributed in a geometrically even 
sense along the trajectory, the angular variable (£ or y) becomes the decision variable in lieu of the flight time. This is 
a non-intuitive variable for the human mission designer to select, however, as any given angular displacement can result 
in dramatically different flight times depending on the geometry of the orbit. Requiring a mission designer to select 
angular displacements instead of flight times is really not feasible. An even more difficult task for a human would 
be to determine, a priori, the temporal spacing between each control node required to achieve perfect geometrical 
distribution of the control nodes. 

Instead of each of these undesirable jobs, in this method we remove the burden of having a mission designer select 
the angular displacement and shift the responsibility to the NLP solver. Then, the universal Kepler solver is used in 
an inverse sense to compute the time-of-flight required to satisfy an even distribution of control points, in addition to 
propagating the state without the need for iteration. The procedure for evaluating a trajectory phase, and for obtaining 
all the necessary partials to compute match point derivatives, is as follows (definitions for the symbols below are 
provided in the Appendix): 


1. Xp is an NLP decision variable and is equal to the total generalized anomaly required to propagate one trajectory 


. . . = xX) 
phase. The size of a propagation segment is Ay; = 4 


2. Starting at the left boundary, propagate inwards to the phase center for + segments of generalized anomaly 
3. For each k" propagation step in the forward half-phase: 


(a) compute a; to determine if the current conic orbit is elliptic, hyperbolic or parabolic (this affects the 
computation of the universal variables U,,, and their derivatives) 


(b) compute the universal functions U,, n =0,1,2,3 


(c) compute sae and 2 Fa n = 0,1, 2,3 (for use in computing step g). 


(d) compute the segment propagation time At; resulting from a propagation through Ax; of generalized 


anomaly as well as oote 


(e) compute the Lagrange coefficients and their derivatives F’,, Fee, kB ko an 7 Gk, 5 et Ren and $+ ah 


(f) compute the propagated state X, = [Te ve] 


and OAtrK 


(g) compute axt |? 


axt ; which form the basis for computing match-point derivatives 


(h) compute the maximum available thrust Tinax,, 


(i) compute Avmax, 
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(j) compute throttle magnitude constraint ||u,|| <1 


(k) apply the impulsive maneuver vi =v, £ugAvimax, 
(1) compute the spacecraft’s mass after the applied impulse mp 
4. Locate the right hand phase boundary using the phase flight time NLP decision variable At, 


5. Repeat step 3 for the backwards half-phase 


6. Se = match point constraint c! and its partials pe (here we denote the phase match point with the dagger 
symbo 


In the procedure above, the - and + superscripts indicate a quantity immediately prior to and after an applied impulse 
respectively. The 3x1 vector uj, contains the control parameters of the k“” segment. The scalar quantity Avmax, 
represents the maximum Av achievable by the spacecraft over the course of the k‘” segment: 


NactiveD Tras: Atr 
Mr 


(1) 


Avmax, = 


In Equation (1), Nactive 1s the number of active thrusters, D is the thruster duty cycle, Tinax, is the maximum available 
thrust for the current segment computed using the propulsion system power model? and my, is the mass of the space- 
craft at the left-hand side of the segment. The spacecraft’s mass across the k*” bounded impulse is computed using 
the following equation: 
J mp1 — |Jug—1|| Armax,_, forward propagation (2) 
— Mrz + |luz]|] Ammax, backward propagation 


where Mmax, = D At 1imax, and At, given by Eq. (11) in the Appendix . Note that the impulsive thrust approxima- 
tion implies mp4 =m, and for notational convenience, we set m;, = Mk. 

One peculiarity of this method is that the maximum size of a maneuver Avmax, can only be determined after 
the propagation time At, has been computed. This is in contrast to using an equivalent time spacing of the nodes 
(MGALT) where segment propagation time is known prior to propagation and hence the trajectory may be computed 
in half-segments. For this reason, maneuvers must occur at the the right-hand boundary of a segment instead of in 
the middle as they would for the original MGALT transcription. This is immaterial, as the location of the impulse is 
more or less arbitrary for an impulsive approximation of continuous-thrust, save for the case of the last segment in any 
half-phase where placing the maneuver on the right-hand boundary results in it occurring exactly at the match-point. 
When both half-phases are considered, this results in the burn at the match-point having a maximum potential size that 
is twice as large as the others. This is not a major concern, however, as the purpose of bounded-impulse methods is to 
produce an accurate estimation of the total mass requirement of a mission, and not necessarily the exact geometry of 
its trajectory. This problem also reduces asymptotically with the use of more control nodes. 

Computation of the match-point constraint partials is carried out in a similar fashion to the original MGALT 
transcription, with the added complexity of the new yp decision variables. 


APPLICATION: RENDEZVOUS WITH COMET 45P/HONDA-MRKOS-PA JDUSAKOVA 


The problem selected to illustrate the benefits of time regularization is a notional rendezvous with the comet 45P/Honda- 
Mrkos-PajduSakova. This is a short-period comet with a period of approximately 5.25 years. We consider a direct flight 
to the comet, with no intermediate flybys. This rendezvous problem was solved once using the MGALT transcription 
(Figure |) and then a second time using the MGALTS transcription (Figure 2) using the exact same configuration 
parameters. The NLP objective function was the minimization of the propellant consumed. 
Specific information about power modeling and throttle logic used is omitted here but is available, and adheres to 
Discovery class mission proposal requirements. !* !4 

Figures | and 2 show that the optimizer arrives at two answers that differ by less than one percent in terms of 
mass delivered to the comet. This is a testament to just how accurate the MGALT transcription is at approximating a 
mission’s mass requirements. The major difference between the two solutions is the flight time required to make the 
transfer. The MGALT trajectory requires 10.10 years, whereas the MGALTS completes the rendezvous in only 7.96 
years. The throttle histories depicted in Figures 3 and 4 shed light on why this is the case. The MGALTS trajectory 
features a much later launch date, and even arrives earlier at the target. It is able to achieve this due to greater number 
of control points located near periapsis for that transcription. The clustering of control nodes there allows the pump- 
up maneuver to happen more efficiently, meaning the spacecraft can match the energy of the comet’s orbit in a much 
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Table 1. 45P/Honda-Mrkos-Pajdusakova rendezvous mission parameters 


Parameter Value 

Maximum allowed initial mass 4000 kg 

Earliest allowed launch date August a7, 2016 

Launch window 5 years 

Maximum flight time unbounded 

Latest allowed arrival date January 1‘, 2030 

Launch vehicle Atlas V (555) NLS-2 

Maximum launch C3 21.7156 km? /s? 

Launch declination bounds [—28.5°, 28.5°] 

Arrival type rendezvous 

Thruster 2 x NEXT TT11 high thrust 

Throttle logic minimum number of thrusters 

Thruster duty cycle 0.95 

Solar power coefficients Yo = 1.32077 71 = —0.10848 
y2 = —0.11665 73 = 0.10843 
ya = —0.01279 

Spacecraft bus power requirement coefficients a, /, = 1.0 
bs/_ = 0.0 
Cs/c = 0.0 

Solar array BOL power (launch at 1 A.U.) 40.0 kW 

Power margin (dpower) 15% 

Post-launch checkout coast 60 days 

Number of segments per phase 200 

Ephemeris SPICE 

SNOPT feasibility tolerance 1.0e-5 

SNOPT optimality tolerance 1.0e-4 

Objective function maximize final mass 


- Event # 1: 

Event # 2: . Launch 

Lt rendezvous Earth 

45P 11/21/2017 
12/29/2027 J C3 = 21.716 km?/s? 
m = 3238 kg DLA = 20.1° 

m = 4000 kg 
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Figure 1. Rendezvous with 45P using MGALT with 200 segments, Atyission = 3690 days 


shorter time. The MGALT solution on the other hand is not able to achieve the same energy change per periapsis 
passage. This works out in a particularly unlucky way for this mission as it results in the spacecraft having to make 
an entire extra revolution of the sun in order to carry out the rendezvous. Since the MGALT node placement is less 
optimal, the optimizer adds an additional revolution in order to place the maneuvers more effectively, which in turn 
requires that the spacecraft launch earlier. 

By increasing the number of control nodes for the MGALT model, a solution in the same family as the MGALTS 
solution can be obtained due to the increased number of nodes at periapsis. Figure 5 shows the new MGALT solution, 
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Event # 1: 


Event # 2: Launch 
Lt rendezvous Earth 


45P : 11/10/2019 
10/28/2027 Cy = 21.716 km?/s? 
m = 3217 kg DLA = 16.6° 

m = 4000 kg 


Figure 2. Rendezvous with 45P using time-regularized MGALTS with 200 segments, Atmission = 2909 days 


— Distance from SUN (AU) 
---- Control magnitude 


Scalar Metric 


01-01-2018 01-01-2019 01-01-2020 01-01-2021 01-01-2022 01-01-2023 01-01-2024 01-01-2025 01-01-2026 01-01-2027 
Epoch 


Figure 3. MGALT throttle magnitude and distance from the sun as a function of time. 


this time using 400 segments instead of 200. Increasing the segment count by degrees will eventually yield a trajectory 
with a more comparable flight time to the MGALTS solution. 
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— _ Distance from SUN (AU) 
---- Control magnitude 


Scalar Metric 


0 pees axa 
01-01-2020 01-01-2021 01-01-2022 01-01-2023 01-01-2024 01-01-2025 01-01-2026 01-01-2027 
Epoch 
Figure 4. MGALTS throttle magnitude and distance from the sun as a function of time. 
1e8 x (km) 
-6 2 2 
T T T T T 
| 2 | Event #1: 
Event # 2: a Launch 
Lt rendezvous Earth 
45P 11/16/2018 
5/19/2028 C3 = 21.716 km?/s” 
m = 3138 kg DLA = 17.4° 
m = 4000 kg 
ioe] 
71-88 
Figure 5. Rendezvous using the time-regularized MGALT transcription with 400 segments 
CONCLUSION 


Time-regularization of the bounded-impulse Kepler low-thrust model has been achieved using an exact, non-iterative 
analytic algorithm. It is also possible to obtain analytic partial derivatives for the nonlinear constraints defining this 
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transcription. This methodology is particularly effective at solving trajectory problems requiring large eccentricity 
changes, especially those requiring multiple revolutions of the central body. By replacing time with the generalized 
anomaly as the independent propagation variable, the control nodes of the MGALT model can be distributed in an 
even geometric fashion, increasing the number of nodes near periapsis. This important feature of the model allows for 
a discrete control history that more closely matches the optimal one that could be found using the necessary conditions 
of the Calculus of Variations. This difference can lead to trajectory solutions with a similar cost function value to 
those produced by the MGALT model, but vary significantly in topology. This is illustrated by the example comet 
rendezvous problem that was solved in this paper where a nearly identical cost function (within 1%) was achieved 
with MGALT and MGALTS, however, the flight time of the MGALTS solution was slightly over two years shorter. 


APPENDIX 
Kepler Propagation with Universal Functions 
The transcription described in this paper utilizes Keplerian two-body propagation using the Lagrange coefficients, | 


re} _ Fr Gr Te-1 
S| 7 E a Rae (3) 


For use with a numerical optimizer, it is generally beneficial to use a two-body propagation method capable of 
robustly propagating any orbit initial condition (i.e. elliptical, parabolic or hyperbolic) initialized by a search method. 
For this reason, a propagation method based on universal orbit variables is employed. The universal variables are 
defined according to the energy regime of the orbit: 


Ui, = xe(1 — yr Sk) 


Us, = 370; 
eso, ee ® (4) 
Us, = X%Sk 
Uo, =1- apU2, 
where, 
1 2 vz 
ae a ee 6) 
Yk = OLN (6) 
1 
Ck = (1 — cos(/#c)) (7) 
1 : 
Sk = a (Vue — sin(/Yr)) (8) 
k 
Uo, = cosh(./—axnXx) 
iy 0 U1, — Jaa sinh(/—aK Xx) (9) 
2h =U _ ox) 
Us, = ac (xe — Ui, ) 
Uo, =1 
U- = 
Gye = 0X 1 (10) 
Uo, = 5U1,,.Xk 
Us, = $U2,.Xk 


The propagation time At; associated with a corresponding change in generalized anomaly x; is computed as follows: 


1 
At, = i (roU1 + a9U2 + Us) (11) 
OAtr 1 ( OU7L OU >) 
= + + 12 
Oxn = ONK ON. OXk oe 
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where, 


rr' Vk 
C= (13) 
VE 
The distance from the central body at the end of the propagation segment is: 
re41 = TRU, + 7KU1, + U2, (14) 
The Lagrange coefficients and their derivatives are given by: 
ree ee (15) 
Tk 
OF ;, 1 OU: 
e= 2s (16) 
OX% Tr OXk 
: Lb 
Ff, = vi Uj, (17) 
TR+ITk 
OF; _ JE Ui, Orns 1 1 OU»i, 
= . (18) 
OXk Tk VRe4d OXk Tk+1 OXK 
1 
G;, = — (r,U, + o,U2 ) (19) 
Ji k k 
dG, 1 ( Oi, +) 
= r +o 20 
OxXk JE * ONE "Oxk a 
G,=1- U2, (21) 
Tk+1 
1 
OG, _ U2. OrK41 OU, (22) 
OXk Trix OXk  Tk+1 OXk 
where, 
0 OU OU: Ou. 
Tk+1 =P OK + op 1x +4 2k (23) 
OXk OXK Oxk — OXk 
Peo = 7RUo, + onU1, + U2, (24) 
Ok4+1 = o~Uo, + (1 = Qk k) Ui, (25) 
OU, On, 
= —azU},; * = Upn_ 1,22 jc. 26 
roma OR aa 1 (26) 
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